Decontam
Decontam
Paths and libraries setting
# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/paths.R")
source("../../../source/functions.R")
# Load supplementary packages
packages <- c("decontam", "kableExtra", "microbiome", "cowplot")
invisible(lapply(packages, require, character.only = TRUE))Detect contaminants
res <- make_decontam(ps, "1D_MED", path_tsv, threshold=0.6)
decontam_asv_MED <- res[[1]]
decontam_df_MED <- res[[2]]
plots_MED <- res[[3]] # 3 contaminants : Pseudomonas, Enhydrobacter, Rahnella1
# print contaminants
decontam_df_MED %>%
kbl() %>%
kable_paper("hover", full_width = F)| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| N0005 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | Enhydrobacter | NA |
| N0315 | Bacteria | Proteobacteria | Gammaproteobacteria | Enterobacterales | Yersiniaceae | Rahnella1 | NA |
| N0852 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Pseudomonadaceae | Pseudomonas | NA |
# histogram of the scores
p <- hist(res[[4]]$p, 100, ylim = c(0,25), xlim = c(0,1), main="", xlab="Score", ylab="Frequency")Save decontam plot
setwd(path_plot)
pdf("1D_MED_decontam_plots.pdf")
for (i in 1:length(plots_MED)) {
print(plots_MED[[i]])
}
dev.off()## quartz_off_screen
## 2
Remove contaminants
# remove contaminants ASV
alltaxa <- taxa_names(ps)
decontam_taxa <- alltaxa[!(alltaxa %in% decontam_asv_MED)]
ps2 <- prune_taxa(decontam_taxa, ps)
# check ps
ps2 <- check_ps(ps2)
ps2## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 67 taxa and 123 samples ]
## sample_data() Sample Data: [ 123 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 67 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 67 reference sequences ]
## [1] 6469784
## [1] 17460
Check blanks
# blanks
ps.blanks <- subset_samples(ps, Location=="Blank")
ps.blanks <- check_ps(ps.blanks)
ps.blanks## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 24 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 24 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 24 reference sequences ]
## [1] 16211
# supprimer blanks
ps.filter <- subset_samples(ps, Location!="Blank")
ps.filter <- check_ps(ps.filter)
ps.filter## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 67 taxa and 115 samples ]
## sample_data() Sample Data: [ 115 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 67 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 67 reference sequences ]
## [1] 6453573
## Compositional = NO2
## 1] Min. number of reads = 952] Max. number of reads = 7361593] Total number of reads = 64535734] Average number of reads = 56118.02608695655] Median number of reads = 251717] Sparsity = 0.6486696950032456] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 16SampleWellPrimer1Primer2LocationFieldCountryOrganSpeciesIndividualIndividualsDateRunControlDnaRead_depth2
## [[1]]
## [1] "1] Min. number of reads = 95"
##
## [[2]]
## [1] "2] Max. number of reads = 736159"
##
## [[3]]
## [1] "3] Total number of reads = 6453573"
##
## [[4]]
## [1] "4] Average number of reads = 56118.0260869565"
##
## [[5]]
## [1] "5] Median number of reads = 25171"
##
## [[6]]
## [1] "7] Sparsity = 0.648669695003245"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
##
## [[8]]
## [1] "8] Number of singletons = 0"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 16"
##
## [[11]]
## [1] "Sample" "Well" "Primer1" "Primer2" "Location"
## [6] "Field" "Country" "Organ" "Species" "Individual"
## [11] "Individuals" "Date" "Run" "Control" "Dna"
## [16] "Read_depth"
Counts
Culex pipiens
## ps_pipiens
ps.pipiens <- subset_samples(ps.filter, Species=="Culex pipiens")
ps.pipiens <- check_ps(ps.pipiens)
ps.pipiens## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 49 taxa and 83 samples ]
## sample_data() Sample Data: [ 83 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 49 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 49 reference sequences ]
## read count: 2574603
ps.pipiens.field <- ps.pipiens %>% subset_samples(Field=="Field")
ps.pipiens.field <- ps.pipiens.field %>% check_ps()
# Culex pipiens - Bosc
cat("ps_pipiens_field_bosc\n\n")## ps_pipiens_field_bosc
ps.pipiens.field.bosc <- ps.pipiens.field %>% subset_samples(Location=="Bosc")
ps.pipiens.field.bosc <- ps.pipiens.field.bosc %>% check_ps()
ps.pipiens.field.bosc## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 36 taxa and 23 samples ]
## sample_data() Sample Data: [ 23 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 36 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 36 reference sequences ]
## read count: 803574
## ps_pipiens_field_camping
ps.pipiens.field.camping <- ps.pipiens.field %>% subset_samples(Location=="Camping Europe")
ps.pipiens.field.camping <- ps.pipiens.field.camping %>% check_ps()
ps.pipiens.field.camping## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 39 taxa and 25 samples ]
## sample_data() Sample Data: [ 25 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 39 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 39 reference sequences ]
## read count: 869617
## ps_pipiens_lavar
ps.pipiens.lavar <- ps.pipiens %>% subset_samples(Location=="Lavar (labo)")
ps.pipiens.lavar <- ps.pipiens.lavar %>% check_ps()
ps.pipiens.lavar## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 48 taxa and 35 samples ]
## sample_data() Sample Data: [ 35 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 48 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 48 reference sequences ]
## read count: 901412
Culex quinquefasciatus
## ps_quinque
ps.quinque <- subset_samples(ps.filter, Species=="Culex quinquefasciatus")
ps.quinque <- check_ps(ps.quinque)
ps.quinque## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 59 taxa and 21 samples ]
## sample_data() Sample Data: [ 21 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 59 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 59 reference sequences ]
## read count: 1925054
## ps_quinque_field
ps.quinque.field <- ps.quinque %>% subset_samples(Field=="Field")
ps.quinque.field <- ps.quinque.field %>% check_ps()
ps.quinque.field## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 47 taxa and 7 samples ]
## sample_data() Sample Data: [ 7 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 47 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 47 reference sequences ]
## read count: 1721053
## ps_quinque_lab
ps.quinque.lab <- ps.quinque %>% subset_samples(Field=="Lab ")
ps.quinque.lab <- ps.quinque.lab %>% check_ps()
ps.quinque.lab## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 45 taxa and 14 samples ]
## sample_data() Sample Data: [ 14 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 45 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 45 reference sequences ]
## read count: 204001
Aedes aegypti
## ps_aedes
ps.aedes <- subset_samples(ps.filter, Species=="Aedes aegypti")
ps.aedes <- check_ps(ps.aedes)
ps.aedes## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 54 taxa and 11 samples ]
## sample_data() Sample Data: [ 11 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 54 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 54 reference sequences ]
## read count: 1953916
Number of reads by sample
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)
metadata_read <- data.frame(ps.filter@sam_data)
make.italic <- function(x) as.expression(lapply(x, function(y) bquote(italic(.(y)))))
levels(metadata_read$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
"Culex pipiens"=make.italic("Culex pipiens"),
"Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))
levels(metadata_read$Location) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", "Wolbachia~'-'~(Slab~TC)")
# metadata_read$f2 <- factor(metadata_read$Species, labels = c("italic(Aedesaegypti)", "italic(Culexpipiens)", "italic(Culexquinquefasciatus)"))
guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 10, face = "italic", colour = "Black", angle = 0)))
p <- ggplot(metadata_read, aes(x = Sample, y = Read_depth))+
geom_bar(position = "dodge", stat = "identity")+
scale_fill_manual(values = col)+
theme_bw() +
theme(axis.text.x = element_text(angle = 90)) +
ggtitle("") +
# guide_italics+
theme(legend.title = element_text(size = 18),
legend.position="bottom",
legend.text=element_text(size=11),
panel.spacing.y=unit(1, "lines"),
panel.spacing.x=unit(0.8, "lines"),
panel.spacing=unit(0,"lines"),
strip.background=element_rect(color="grey30", fill="grey90"),
panel.border=element_rect(color="grey90"),
axis.ticks.x=element_blank(),
strip.text.x = element_text(size = 12)) +
facet_wrap(~Species+Location+Organ, scales = "free_x", ncol=4, labeller=label_parsed)+
labs(y="Sequence counts")+
ylim(0, 550000)+
geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=2, hjust=-0.25, vjust=0.25, angle=90)+
guides(fill=guide_legend(title="Oligotype", label.theme = element_text(size = 10, face = "italic", colour = "Black", angle = 0)))
pRemove sample with number of reads <1000
a <- prune_samples(sample_sums(ps.filter)<=1000, ps.filter)
test <- data.frame(a@sam_data)
test <- test[test$Sample!="NP20" & test$Sample!="NP29" & test$Sample!="NP30" & test$Sample!="NP34" & test$Sample!="NP36",]
toremove <- test$Sample
ps.filter <- subset_samples(ps.filter, !(Sample %in% toremove))
ps.filter## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 67 taxa and 113 samples ]
## sample_data() Sample Data: [ 113 samples by 16 sample variables ]
## tax_table() Taxonomy Table: [ 67 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 67 reference sequences ]
## [1] 6452623
Final number of reads by sample
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)
metadata_read <- data.frame(ps.filter@sam_data)
levels(metadata_read$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
"Culex pipiens"=make.italic("Culex pipiens"),
"Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))
levels(metadata_read$Location) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"),"- (Slab TC)")))
metadata_read_whole <- metadata_read[metadata_read$Organ=="Whole",]
metadata_read_ovary <- metadata_read[metadata_read$Organ=="Ovary",]
guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 10, face = "italic", colour = "Black", angle = 0)))
p <- ggplot(metadata_read_whole, aes(x = Sample, y = Read_depth))+
geom_bar(position = "dodge", stat = "identity")+
scale_fill_manual(values = col)+
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust=1, hjust=1, size=15)) +
ggtitle("") +
guide_italics+
theme(legend.title = element_text(size = 18),
legend.position="bottom",
legend.text=element_text(size=16),
panel.spacing.y=unit(1, "lines"),
panel.spacing.x=unit(0.8, "lines"),
panel.spacing=unit(0,"lines"),
strip.background=element_rect(color="grey30", fill="grey90"),
panel.border=element_rect(color="grey90"),
axis.ticks.x=element_blank(),
strip.text.x = element_text(size = 16),
axis.text = element_text(size = 15)) +
facet_wrap(~Species+Location+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
labs(y="Sequence counts")+
ylim(0, 1000000)+
#geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=4, hjust=-0.25, vjust=0.25, angle=90)+
geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=4, angle=90, hjust=-0.1, vjust=0.25)+
guides(fill=guide_legend(title="Oligotype"))
pp2 <- ggplot(metadata_read_ovary, aes(x = Sample, y = Read_depth))+
geom_bar(position = "dodge", stat = "identity")+
scale_fill_manual(values = col)+
theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust=1, hjust=1, size=15)) +
ggtitle("") +
guide_italics+
theme(legend.title = element_text(size = 18),
legend.position="bottom",
legend.text=element_text(size=16),
panel.spacing.y=unit(1, "lines"),
panel.spacing.x=unit(0.8, "lines"),
panel.spacing=unit(0,"lines"),
strip.background=element_rect(color="grey30", fill="grey90"),
panel.border=element_rect(color="grey90"),
axis.ticks.x=element_blank(),
strip.text.x = element_text(size = 16),
axis.text = element_text(size = 15)) +
facet_wrap(~Species+Location+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
labs(y="Sequence counts")+
ylim(0, 1000000)+
geom_text(aes(label=Read_depth), position=position_dodge(width=0.5), size=4, angle=90, hjust=-0.1, vjust=0.25)+
guides(fill=guide_legend(title="Oligotype"))
p2# panels
p_group <- plot_grid(p+theme(legend.position="none"),
p2+theme(legend.position="none"),
nrow=2,
ncol=1)+
draw_plot_label(c("A", "B"), c(0, 0), c(1, 0.5), size = 15)
legend_plot <- get_legend(p + theme(legend.position="bottom"))
p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_countsFinal number of reads by Genus
new_names_genus <- c("Wolbachia",
"Asaia",
"Legionella",
"Elizabethkingia",
"Chryseobacterium",
"Erwinia",
"Morganella",
"Pseudomonas",
"Delftia",
"Methylobacterium-Methylorubrum",
"Serratia",
"Coetzeea",
"NA"
)
col_genus <- c("Wolbachia"="#FEB24C",
"Asaia"="#10E015",
"Legionella"="#DE3F23",
"Elizabethkingia"="#66A7ED",
"Chryseobacterium"="#F899FF",
"Erwinia"="#FFE352",
"Morganella"="#F5E4D3",
"Pseudomonas"="#DBF5F0",
"Delftia"="#C7C5B7",
"Methylobacterium-Methylorubrum"="blue",
"Serratia"="#B136F5",
"Coetzeea"="red",
"NA"="grey")
p <- plot_bar(ps.filter, "Genus", "Abundance", "Genus")+
geom_bar(aes(color=Genus, fill=Genus), stat="identity", position="stack")+
scale_fill_manual(values = col_genus)+
labs(x="Genus", y="Number of read") +
facet_wrap(~sample_Species+Location+Organ, scales = "free", ncol=3)## Warning in psmelt(physeq): The sample variables:
## Sample
## have been renamed to:
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
## Warning in psmelt(physeq): The sample variables:
## Species
## have been renamed to:
## sample_Species
## to avoid conflicts with taxonomic rank names.
data <- p$data
labels = c("Wolbachia"=make.italic("Wolbachia"),
"Asaia"=make.italic("Asaia"),
"Legionella"=make.italic("Legionella"),
"Elizabethkingia"=make.italic("Elizabethkingia"),
"Chryseobacterium"=make.italic("Chryseobacterium"),
"Erwinia"=make.italic("Erwinia"),
"Morganella"=make.italic("Morganella"),
"Pseudomonas"=make.italic("Pseudomonas"),
"Delftia"=make.italic("Delftia"),
"Methylobacterium-Methylorubrum"=make.italic("Methylobacterium-Methylorubrum"),
"Serratia"=make.italic("Serratia"),
"Coetzeea"=make.italic("Coetzeea"),
"NA"
)
data_pipiens_whole <- data[data$sample_Species=="Culex pipiens" & data$Organ=="Whole",]
data_pipiens_ovary <- data[data$sample_Species=="Culex pipiens" & data$Organ=="Ovary",]
data_quinque_whole <- data[data$sample_Species=="Culex quinquefasciatus" & data$Organ=="Whole",]
data_quinque_ovary <- data[data$sample_Species=="Culex quinquefasciatus" & data$Organ=="Ovary",]
data_aedes_whole <- data[data$sample_Species=="Aedes aegypti" & data$Organ=="Whole",]
data_aedes_ovary <- data[data$sample_Species=="Aedes aegypti" & data$Organ=="Ovary",]
# pipiens
## whole
### all
pipiens1 <- percent_taxonomic_plot_test(data=data_pipiens_whole, variable=sample_Species, selection="Culex pipiens", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(W)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'sample_Species' (override with `.groups` argument)
### field
pipiens2 <- percent_taxonomic_plot_test(data=data_pipiens_whole, variable=Location, selection="Bosc", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(W)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'Location' (override with `.groups` argument)
## Warning: Unknown or uninitialised column: `sample_Species`.
## Warning: Unknown or uninitialised column: `Organ`.
pipiens3 <- percent_taxonomic_plot_test(data=data_pipiens_whole, variable=Location, selection="Camping Europe", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(W)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'Location' (override with `.groups` argument)
## Warning: Unknown or uninitialised column: `sample_Species`.
## Warning: Unknown or uninitialised column: `Organ`.
### labo
pipiens4 <- percent_taxonomic_plot_test(data=data_pipiens_whole, variable=Location, selection="Lavar (labo)", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(W)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'Location' (override with `.groups` argument)
## Warning: Unknown or uninitialised column: `sample_Species`.
## Warning: Unknown or uninitialised column: `Organ`.
## ovary
## all
pipiens5 <- percent_taxonomic_plot_test(data=data_pipiens_ovary, variable=sample_Species, selection="Culex pipiens", group="Genus", organ="Ovary", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(O)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'sample_Species' (override with `.groups` argument)
# quinque
## whole
### all
quinque1 <- percent_taxonomic_plot_test(data=data_quinque_whole, variable=sample_Species, selection="Culex quinquefasciatus", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(W)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'sample_Species' (override with `.groups` argument)
### field
quinque2 <- percent_taxonomic_plot_test(data=data_quinque_whole, variable=Location, selection="Guadeloupe", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(W)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'Location' (override with `.groups` argument)
## Warning: Unknown or uninitialised column: `sample_Species`.
## Warning: Unknown or uninitialised column: `Organ`.
### labo
quinque3 <- percent_taxonomic_plot_test(data=data_quinque_whole, variable=Location, selection="Wolbachia -",group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(W)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'Location' (override with `.groups` argument)
## Warning: Unknown or uninitialised column: `sample_Species`.
## Warning: Unknown or uninitialised column: `Organ`.
quinque3[["labels"]][["title"]] <- expression(paste(italic("Wolbachia"), "- (Slab TC)"))
## ovary
## all
quinque4 <- percent_taxonomic_plot_test(data=data_quinque_ovary, variable=sample_Species, selection="Culex quinquefasciatus",group="Genus", organ="Ovary", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(O)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'sample_Species' (override with `.groups` argument)
# aedes
## whole
### all
aedes1 <- percent_taxonomic_plot_test(data=data_aedes_whole, variable=sample_Species, selection="Aedes aegypti",group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(W)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'sample_Species' (override with `.groups` argument)
## ovary
## all
aedes2 <- percent_taxonomic_plot_test(data=data_aedes_ovary, variable=sample_Species, selection="Aedes aegypti", group="Genus", organ="Ovary", new_names=new_names_genus, col=col_genus, mylabels=labels)+
labs(tag = "(O)")+
theme(plot.tag.position = "topright")## `summarise()` regrouping output by 'sample_Species' (override with `.groups` argument)
p_genus <- plot_grid(pipiens1+ theme(legend.position="none"),
pipiens2+ theme(legend.position="none"),
pipiens3+ theme(legend.position="none"),
pipiens4+ theme(legend.position="none"),
pipiens5+ theme(legend.position="none"),
quinque1+ theme(legend.position="none"),
quinque2+ theme(legend.position="none"),
quinque3+ theme(legend.position="none"),
quinque4+ theme(legend.position="none"),
plot.new(),
aedes1+ theme(legend.position="none"),
aedes2+ theme(legend.position="none"),
plot.new(),
plot.new(),
plot.new(),
nrow=3,
ncol=5)+
draw_plot_label(c("A", "B", "C"), c(0, 0, 0), c(1, 2/3, 1/3), size = 15)Counts after removing samples
Create a dataframe
x <- c("Culex pipiens", "Bosc", "Camping Europe", "Lavar (lab)",
"Culex quinquefasciatus", "AGuadeloupe", "Wolbachia - (Slab TC)",
"Aedes aegypti", "BGuadeloupe",
"Total")
y <- c("Reads", "Oligotypes", "Samples")
df <- data.frame(matrix(ncol=3, nrow=10))
rownames(df) <- x
colnames(df) <- y
df2 <- df
df3 <- dfWhole + Ovary
# Culex pipiens
ps.pipiens <- subset_samples(ps.filter, Species=="Culex pipiens")
ps.pipiens <- check_ps(ps.pipiens)
## All location
### oligotype
nrow(ps.pipiens@otu_table) -> df["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data) -> df["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens) -> df["Culex pipiens", "Reads"]
## Bosc
### oligotype
nrow(otu_table(ps.pipiens %>%
subset_samples(Location=="Bosc") %>%
check_ps())) -> df["Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Location=="Bosc")) -> df["Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Location=="Bosc")) -> df["Bosc", "Reads"]
## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens %>%
subset_samples(Location=="Camping Europe") %>%
check_ps())) -> df["Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Location=="Camping Europe")) -> df["Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Location=="Camping Europe")) -> df["Camping Europe", "Reads"]
## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens %>%
subset_samples(Location=="Lavar (labo)") %>%
check_ps())) -> df["Lavar (lab)", "Oligotypes"]
### samples
nrow(ps.pipiens@sam_data %>% subset_samples(Location=="Lavar (labo)")) -> df["Lavar (lab)", "Samples"]
### reads
compute_read_counts(ps.pipiens %>% subset_samples(Location=="Lavar (labo)")) -> df["Lavar (lab)", "Reads"]
# Culex quinquefasciatus
ps.quinque <- subset_samples(ps.filter, Species=="Culex quinquefasciatus")
ps.quinque <- check_ps(ps.quinque)
## All location
### oligotype
nrow(ps.quinque@otu_table) -> df["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data) -> df["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque) -> df["Culex quinquefasciatus", "Reads"]
## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque %>%
subset_samples(Location=="Guadeloupe") %>%
check_ps())) -> df["AGuadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Location=="Guadeloupe")) -> df["AGuadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Location=="Guadeloupe")) -> df["AGuadeloupe", "Reads"]
## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque %>%
subset_samples(Location=="Wolbachia -") %>%
check_ps())) -> df["Wolbachia - (Slab TC)", "Oligotypes"]
### samples
nrow(ps.quinque@sam_data %>% subset_samples(Location=="Wolbachia -")) -> df["Wolbachia - (Slab TC)", "Samples"]
### reads
compute_read_counts(ps.quinque %>% subset_samples(Location=="Wolbachia -")) -> df["Wolbachia - (Slab TC)", "Reads"]
# Aedes aegypti
ps.aedes <- subset_samples(ps.filter, Species=="Aedes aegypti")
ps.aedes <- check_ps(ps.aedes)
## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes %>%
subset_samples(Location=="Guadeloupe") %>%
check_ps())) -> df[c("BGuadeloupe","Aedes aegypti"), "Oligotypes"]
### samples
nrow(ps.aedes@sam_data %>% subset_samples(Location=="Guadeloupe")) -> df[c("BGuadeloupe","Aedes aegypti"), "Samples"]
### reads
compute_read_counts(ps.aedes %>% subset_samples(Location=="Guadeloupe")) -> df[c("BGuadeloupe","Aedes aegypti"), "Reads"]
# Total
### oligotype
nrow(ps.filter@otu_table) -> df["Total", "Oligotypes"]
### samples
nrow(ps.filter@sam_data) -> df["Total", "Samples"]
### reads
compute_read_counts(ps.filter) -> df["Total", "Reads"]Whole
# Culex pipiens
ps.pipiens.whole <- ps.pipiens %>% subset_samples(Organ=="Whole")
ps.pipiens.whole <- ps.pipiens.whole %>% check_ps()
## All location
### oligotype
nrow(ps.pipiens.whole@otu_table) -> df2["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data) -> df2["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole) -> df2["Culex pipiens", "Reads"]
## Bosc
### oligotype
nrow(otu_table(ps.pipiens.whole %>%
subset_samples(Location=="Bosc") %>%
check_ps())) -> df2["Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Location=="Bosc")) -> df2["Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Location=="Bosc")) -> df2["Bosc", "Reads"]
## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.whole %>%
subset_samples(Location=="Camping Europe") %>%
check_ps())) -> df2["Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Location=="Camping Europe")) -> df2["Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Location=="Camping Europe")) -> df2["Camping Europe", "Reads"]
## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.whole %>%
subset_samples(Location=="Lavar (labo)") %>%
check_ps())) -> df2["Lavar (lab)", "Oligotypes"]
### samples
nrow(ps.pipiens.whole@sam_data %>% subset_samples(Location=="Lavar (labo)")) -> df2["Lavar (lab)", "Samples"]
### reads
compute_read_counts(ps.pipiens.whole %>% subset_samples(Location=="Lavar (labo)")) -> df2["Lavar (lab)", "Reads"]
# Culex quinquefasciatus
ps.quinque.whole <- subset_samples(ps.quinque, Organ=="Whole")
ps.quinque.whole <- check_ps(ps.quinque.whole)
## All location
### oligotype
nrow(ps.quinque.whole@otu_table) -> df2["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data) -> df2["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.whole) -> df2["Culex quinquefasciatus", "Reads"]
## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.whole %>%
subset_samples(Location=="Guadeloupe") %>%
check_ps())) -> df2["AGuadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Location=="Guadeloupe")) -> df2["AGuadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Location=="Guadeloupe")) -> df2["AGuadeloupe", "Reads"]
## Wolbachia -
### oligotype
nrow(otu_table(ps.quinque.whole %>%
subset_samples(Location=="Wolbachia -") %>%
check_ps())) -> df2["Wolbachia - (Slab TC)", "Oligotypes"]
### samples
nrow(ps.quinque.whole@sam_data %>% subset_samples(Location=="Wolbachia -")) -> df2["Wolbachia - (Slab TC)", "Samples"]
### reads
compute_read_counts(ps.quinque.whole %>% subset_samples(Location=="Wolbachia -")) -> df2["Wolbachia - (Slab TC)", "Reads"]
# Aedes aegypti
ps.aedes.whole <- subset_samples(ps.aedes, Organ=="Whole")
ps.aedes.whole <- check_ps(ps.aedes.whole)
## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.whole %>%
subset_samples(Location=="Guadeloupe") %>%
check_ps())) -> df2[c("BGuadeloupe","Aedes aegypti"), "Oligotypes"]
### samples
nrow(ps.aedes.whole@sam_data %>% subset_samples(Location=="Guadeloupe")) -> df2[c("BGuadeloupe","Aedes aegypti"), "Samples"]
### reads
compute_read_counts(ps.aedes.whole %>% subset_samples(Location=="Guadeloupe")) -> df2[c("BGuadeloupe","Aedes aegypti"), "Reads"]
# Total
ps.whole <- ps.filter %>% subset_samples(Organ=="Whole")
ps.whole <- ps.whole %>% check_ps()
### oligotype
nrow(ps.whole@otu_table) -> df2["Total", "Oligotypes"]
### samples
nrow(ps.whole@sam_data) -> df2["Total", "Samples"]
### reads
compute_read_counts(ps.whole) -> df2["Total", "Reads"]Ovary
# Culex pipiens
ps.pipiens.ovary <- ps.pipiens %>% subset_samples(Organ=="Ovary")
ps.pipiens.ovary <- ps.pipiens.ovary %>% check_ps()
## All location
### oligotype
nrow(ps.pipiens.ovary@otu_table) -> df3["Culex pipiens", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data) -> df3["Culex pipiens", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary) -> df3["Culex pipiens", "Reads"]
## Bosc
### oligotype
nrow(otu_table(ps.pipiens.ovary %>%
subset_samples(Location=="Bosc") %>%
check_ps())) -> df3["Bosc", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Location=="Bosc")) -> df3["Bosc", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Location=="Bosc")) -> df3["Bosc", "Reads"]
## Camping Europe
### oligotype
nrow(otu_table(ps.pipiens.ovary %>%
subset_samples(Location=="Camping Europe") %>%
check_ps())) -> df3["Camping Europe", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Location=="Camping Europe")) -> df3["Camping Europe", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Location=="Camping Europe")) -> df3["Camping Europe", "Reads"]
## Lavar (labo)
### oligotype
nrow(otu_table(ps.pipiens.ovary %>%
subset_samples(Location=="Lavar (labo)") %>%
check_ps())) -> df3["Lavar (lab)", "Oligotypes"]
### samples
nrow(ps.pipiens.ovary@sam_data %>% subset_samples(Location=="Lavar (labo)")) -> df3["Lavar (lab)", "Samples"]
### reads
compute_read_counts(ps.pipiens.ovary %>% subset_samples(Location=="Lavar (labo)")) -> df3["Lavar (lab)", "Reads"]
# Culex quinquefasciatus
ps.quinque.ovary <- subset_samples(ps.quinque, Organ=="Ovary")
ps.quinque.ovary <- check_ps(ps.quinque.ovary)
## All location
### oligotype
nrow(ps.quinque.ovary@otu_table) -> df3["Culex quinquefasciatus", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data) -> df3["Culex quinquefasciatus", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary) -> df3["Culex quinquefasciatus", "Reads"]
## Guadeloupe
### oligotype
nrow(otu_table(ps.quinque.ovary %>%
subset_samples(Location=="Guadeloupe") %>%
check_ps())) -> df3["AGuadeloupe", "Oligotypes"]
### samples
nrow(ps.quinque.ovary@sam_data %>% subset_samples(Location=="Guadeloupe")) -> df3["AGuadeloupe", "Samples"]
### reads
compute_read_counts(ps.quinque.ovary %>% subset_samples(Location=="Guadeloupe")) -> df3["AGuadeloupe", "Reads"]
## Wolbachia -
### oligotype
# nrow(otu_table(ps.quinque.ovary %>%
# subset_samples(Location=="Wolbachia -") %>%
# check_ps())) -> df3["Wolbachia - (Slab TC)", "Oligotypes"]
### samples
# nrow(ps.quinque.ovary@sam_data %>% subset_samples(Location=="Wolbachia -")) -> df3["Wolbachia - (Slab TC)", "Samples"]
# ### reads
# compute_read_counts(ps.quinque.ovary %>% subset_samples(Location=="Wolbachia -")) -> df3["Wolbachia - (Slab TC)", "Reads"]
# Aedes aegypti
ps.aedes.ovary <- subset_samples(ps.aedes, Organ=="Ovary")
ps.aedes.ovary <- check_ps(ps.aedes.ovary)
## Guadeloupe
### oligotype
nrow(otu_table(ps.aedes.ovary %>%
subset_samples(Location=="Guadeloupe") %>%
check_ps())) -> df3[c("BGuadeloupe","Aedes aegypti"), "Oligotypes"]
### samples
nrow(ps.aedes.ovary@sam_data %>% subset_samples(Location=="Guadeloupe")) -> df3[c("BGuadeloupe","Aedes aegypti"), "Samples"]
### reads
compute_read_counts(ps.aedes.ovary %>% subset_samples(Location=="Guadeloupe")) -> df3[c("BGuadeloupe","Aedes aegypti"), "Reads"]
# Total
ps.ovary <- ps.filter %>% subset_samples(Organ=="Ovary")
ps.ovary <- ps.ovary %>% check_ps()
### oligotype
nrow(ps.ovary@otu_table) -> df3["Total", "Oligotypes"]
### samples
nrow(ps.ovary@sam_data) -> df3["Total", "Samples"]
### reads
compute_read_counts(ps.ovary) -> df3["Total", "Reads"]
df3[is.na(df3)] <- 0